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While the rheology of non-Brownian suspensions in the dilute regime 
is well-understood, their behavior in the dense limit remains mysti- 
fying. As the packing fraction of particles increases, particle motion 
becomes more collective, leading to a growing length scale and scal- 
ing properties in the rheology as the material approaches the jam- 
ming transition. There is no accepted microscopic description of this 
phenomenon. However, in recent years it has been understood that 
the elasticity of simple amorphous solids is governed by a critical 
point, the unjamming transition where the pressure vanishes, and 
where elastic properties display scaling and a diverging length scale. 
The correspondence between these two transitions is at present un- 
clear. Here we show that for a simple model of dense flow, which 
we argue captures the essential physics near the jamming threshold, 
a formal analogy can be made between the rheology of the flow 
and the elasticity of simple networks. This analogy leads to a new 
conceptual framework to relate microscopic structure to rheology. It 
enables us to define and compute numerically normal modes and a 
density of states. We find striking similarities between the density of 
states in flow, and that of amorphous solids near unjamming: both 
display a plateau above some frequency scale ui* ~ \z c — z\, where 
z is the coordination of the network of particles in contact, z c = 2D 
where D is the spatial dimension. However, a spectacular difference 
appears: the density of states in flow displays a single mode at an- 
other frequency scale w m i n << lu* governing the divergence of the 
viscosity. 
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Introduction 

Suspensions are heterogeneous fluids containing solid parti- 
cles. Their viscosity in the dilute regime was computed early 
on by Einstein and Batchelor [1]. However, as the packing 
fraction cf> increases, steric hindrance becomes dominant and 
particles move under stress in an increasingly coordinated way. 
For non-Brownian particles, the viscosity diverges as the sus- 
pension jams into an amorphous solid. This jamming tran- 
sition is a non-equilibrium critical phenomena: the rheology 
displays scaling laws [2, 3, 4, 5, 6, 7, 8] and a growing length 
scale [2, 3, 9, 10]. Jamming occurs more generally in driven 
materials made of repulsive particles, such as in aerial granu- 
lar flows [11] with similar rheological phenomenology [12], and 
where large eddies are observed as jamming is approached [9]. 
This phenomenology bears similarity with that of the glass 
transition, where steric hindrance increases upon cooling, and 
where the dynamics becomes increasingly collective as relax- 
ation times grow [13]. In dense granular materials and in 
supercooled liquids, the physical origin of collective dynamics 
and associated rheological phenomena remain elusive. 

Recent progress has been made on a related problem, the 
unjamming transition where a solid made of repulsive soft par- 
ticles is isotropically decompressed toward vanishing pressure 
[14, 15, 16, 17, 18, 19, 20]. In this situation various properties 
of the amorphous solid, such as elasticity, transport and force 
propagation, display scaling with the distance from thresh- 
old [14, 21, 22] and are characterized by a diverging length 
scale [17, 18]. The packing geometry also displays scaling: 
the average number of contacts or coordination z of jammed 
particles, is observed to converge [5] to the minimal value al- 



lowing mechanical stability z c derived by Maxwell [23]. The- 
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Fig. 1. (a) Snapshot of a configuration visited under flow from simulations in two 
dimensions (see Methods Section for details). Small particles are in blue, large parti- 
cles are in yellow. The width of the black segments connecting particles is proportional 
to the contact force. The arrows indicate the shear direction, (b) Example of the 
evolution of the stress with time. Instantaneous jumps upward correspond to the cre- 
ation of new contacts. Stress relaxes smoothly however within periods where contacts 
do not change. Interestingly, this situation is opposite to the plasticity of amorphous 
elastic solids, where stress loads continuously before relaxing by sudden plastic events, 
as illustrated in (c) for compressed soft elastic particles under quasi-static shear. 
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oretically, these observations can be unified by the realiza- 
tion that the vibrational spectrum of these simple amorphous 
solids must display a cross-over at some threshold frequency 
[15], distinct from Anderson localization, above which vibra- 
tional modes are still extended but poorly transport energy 
[16]. This threshold frequency, corresponding to the so-called 
boson peak ubiquitously observed in glasses [24], is governed 
by coordination and applied pressure [25] and vanishes at the 
unjamming threshold. 

Under an imposed shear flow, materials near the jamming 
transition are strongly anisotropic and particles receive a net 
force from particles they are in contact with. When an amor- 
phous solid is decompressed toward its unjamming transition 
however, configurations are isotropic and forces are balanced. 
Despite these significant differences, both transitions deal with 
the emergence or disappearance of rigidity, and it is natural 
to seek a common description of these phenomena. Such a 
framework would unify the elastic and vibrational properties 
of amorphous solids, much studied in the context of granu- 
lar materials, glasses [24] and gels of semi-flexible polymers 
[26], with the rheology of suspension flows. Herein we build 
the foundation of this framework, by focusing on a model of 
a suspension flow of non-Brownian hard particles where hy- 
drodynamic interactions are neglected [2]. The relationship 
between this model and real suspensions is not obvious. Here 
we furnish numerical observations and physical arguments to 
support the contention that this model captures the relevant 




10° 10 1 10 z 10 3 10 4 10 s 



p 

Fig. 2. The relation (ftip) within ASM is computed by averaging the normalized 
pressure p at different packing fractions. This relation is well-captured by Eq.(l) 
introduced to describe real suspension flows, with (a) <p c S3 0.842 and a 0.46 
for D = 2 and N = 4096 particles, (b) <j) c « 0.647 and a S3 0.38 for D = 3 
and JV = 1000. The friction law fj,(p) follows Eq.(2), with (a) K, 0.093 and 
(8 « 0.35 for D = 2, and (b) fl c RJ 0.115 and /3 « 0.33 for D = 3. The 
insets display — [l c averaged at fixed pressure. The main plots show non-averaged 
data, and colors label the different packing fractions were measurements were made, 
growing from blue to red (see Methods Section). 



physics near jamming. The next step is to show that for this 
model, we can build a formal analogy between rheology and 
elasticity. This analogy allows us to simulate accurately and 
to measure with high precision key aspects of the particle orga- 
nization, in particular the packing coordination. More funda- 
mentally, it enables us to define and compute, in flows of hard 
particles, normal modes and their density, which are the nat- 
ural mathematical objects connecting geometry to rheology. 
We find striking similarities between the jamming and the un- 
jamming transitions: the coordination converges to Maxwell's 
value with a non-trivial scaling with pressure. Furthermore, 
the density of states displays a transition that is characteris- 
tic of vibrational properties of amorphous solids. However, a 
fundamental difference appears: the density of states displays 
bi-scaling instead of the simple scaling found near unjamming, 
with one mode at low- frequency that strongly couples to shear, 
which is responsible for the rapid divergence in viscosity. 

Constitutive Relations 

Dimensional analysis implies that suspensions of over-damped 
hard particles are controlled by one dimensionless parameter 
[6], the viscous number I = 7^0 /p where p is the particle 
pressure, 770 is the fluid viscosity and 7 is the strain rate. 
For convenience we shall consider its inverse, the normalized 
pressure p = p/'yijo. The steady-state rheology is then en- 
tirely determined by two constitutive relations that allow the 
computation of the flow in various geometries [6, 12]: 4>(p) 
describes the dilatency of the material, and = a/p its 

effective friction, where a is the shear stress. Near jamming 
it is observed experimentally that: 

0c-<Mp)«l/p Q , [1] 
fi(p) - fj, c oc 1//, [2] 

with a w 1/2 [6] and w 0.4 [10, 6, 27]. Currently these 
constitutive relations are phenomenological, and lack a mi- 
croscopic theory. 

The Affine Solvent Model 

We consider a model of hard frictionless particles immersed in 
a viscous fluid, that we call the Affine Solvent Model (ASM). 
The dynamics is over-damped, and hydrodynamic interactions 
are neglected: the viscous drag on a particle is equal to the 
difference between the imposed velocity of the underlying fluid 
and the particle velocity, times a coefficient £o- The flow of 
the fluid phase is undisturbed by the particles and is chosen 
to be an affine shear of strain rate 7. For hard particles at 
fixed packing fraction, changing 7 simply rescales time and 
the stress tensor must therefore be proportional to 7 [30]. 
Note that for this model we can renormalize the pressure as 
p = p/7^0 in two dimensions and p = pd/'y^o'm three dimen- 
sions, where d is taken to be the diameter of the small particles 
in our bi-disperse numerics. 

We constructed an event-driven simulation scheme (see 
[31] and the Methods Section below) based on an exact equa- 
tion of motion derived below, and on the actualization of the 
contact network, to simulate flow while carefully extracting 
the geometry of configurations visited. In Fig.(l) we present 
a snapshot of a configuration from our simulations in two- 
dimensions, and a typical stress-strain signal. A movie (Video 
SI) of our simulations can be found in the Supplementary In- 
formation. 

Our results for the constitutive relations are shown in 
Fig. (2). We first note that fluctuations in pressure for a 
given (/>, proportional to the extent along the p-axis of the 
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corresponding data-points cloud, diverge as jamming is ap- 
proached, as expected from the observation that 4> c fluctuates 
in a finite size system [35]. This finite size effect creates a 
difficulty to accurately extract the exponent a that charac- 
terizes the dilatancy, see Eq.(l). Our simulation method does 
not allow for the study of larger systems, which would im- 
prove our estimation of this exponent; nevertheless, our re- 
sults in two-dimensions are in good agreement with previous 
work [2, 28, 29]. Our results for three dimensions are similar 
to experimental observations: in three dimensions the friction 
law follows Eq.(2) with « 0.33, whereas dilatency follows 
Eq.(l) with a « 0.38. These exponents are close to the ex- 
perimentally observed exponents, supporting the claim that 
hydrodynamical interactions have weak effects, if any, on the 
critical behavior near jamming. Note that the exponents are 
found to be similar in two and three dimensions, raising the 
possibility that they are actually the same. 

In sharp contrast with the fluctuations of pressure for a 
fixed packing fraction, the fluctuations in various quantities 
considered as a function of pressure are very limited, and 
seem to remain regular throughout the entire range of pres- 
sures simulated, see Figs. (2), (3) and (5). In addition, we 
detect no finite size effects when computing means as a func- 
tion of pressure, see, for example, Fig. (5). Thus pressure is a 
much more suitable variable than packing fraction to quantify 
the distance from threshold and to study criticality in a finite 
system. 



Relationship between dynamics and packing geometry 

Within ASM the viscous drag force Fk acting on particle k at 
position Rt is Ft = — £o(Vk — v(Rk)), where Vk and V c (Rk) 
are the particle and fluid velocities respectively, and £o char- 
acterizes the viscous drag. This linear relation can be written 
for all particles in compact notation: 

\F) = -Zo(\V)-\V l )) = -£o\Vn, [3] 
where the upper case ket notation \X) indicates a vector of di- 
mension ND, D being the spatial dimension and N the num- 
ber of particles. |V na ) = \V) — \ V f ) is the so-called non-affine 
velocity. 

In addition to the viscous drag, a force fikHik is ex- 
erted on the k th particle from all particles i which are in 
contact with particle k, where Hik is the unit vector along 
Rk — Ri- In the viscous limit considered here, forces are bal- 
anced F k + Yl fikftik — 0, where the sum is over all particles 
i in contact with k. This equation can be written in compact 
notation: 

|^>+T|/>=0, [4] 

where the lower case ket notation \x) indicates a vector of di- 
mension iV c , the number of contacts. The linear operator T 
is of dimension ND x N c , its non-zero elements correspond to 
the vectors Hik, for particles i and k in contact. 

The last condition defining ASM is that particles cannot 
overlap. This condition is best expressed by considering the 
network of contacts, as illustrated in Fig.(l). This network 
rewires by instantaneous events where contacts open or close 
(see Video SI in the S.I). In between these discrete events, 
the network is conserved, and the distance between particles 
in contact is fixed. This condition implies that 

(V k - Vi) ■ n lk = , [5] 

for all contacts. The N c linear constraints in the form of 
Eq.(5) (one constraint for each pair of particles in contact) 
can be written in compact notation as: 

S\V) = S\V f ) + S\V™} = or S\V na ) = S\V f ), [6] 



where 5 is the N c x ND linear operator that computes the 
change in pairwise distances \6r) between particles in contact, 
induced by a displacement of the particles \SR): S\SR) = \6r). 
A direct inspection of the elements of 5 indicates that it is 
the transpose of T: 5 = T* . 

Now we may derive an expression for the viscosity. Oper- 
ating with 5 on Eq.(3), and using Eqs.(4) and (6), we obtain: 



Af\f) = -£ 5|V f > 



[7] 



where Af = ST = 55* is a N c x iV c symmetric operator, which 
is generally invertible when the viscosity is finite (see below). 
Af contains the information on the topology of the contact 
network (who is in contact with whom) and the contact ori- 
entations. Af characterizes how a contact force field |/) is un- 
balanced {f\N\f) = (/|T*T|/> = (F\F). The term 5|V f ) on 
the right hand-side of Eq.(7) is not singular near jamming, as 
it is the rate at which the contacts lengths would change (gen- 
erating gaps or overlaps between particles) if particles were 
following the affine flow of the fluid. One can thus introduce 
the notation S\V f ) = 7 | 7 ). For a pure shear in the (x,y) 
plane, the components of I7) are 7^ = rij(Hij ■ e x )(ftij ■ e y ), 
where rij is the distance between particles i and j. Inverting 
Eq.(7) and using this notation we derive an expression for the 
contact forces: 

|/) = -^07^- 1 |7> [8] 

This result, together with Eqs.(3) and (4) yields an expression 
for the dynamics: |V na ) = — A /S t Af~ 1 17). The shear stress a 
carried by the particles is related to the contact forces by the 
relation [36] a = -e x -{Yjij fijnjHij igiHij) ■ e y /fl = -{f\^f)/Q, 
where f2 is the volume of the system. We thus obtain from 
Eq.(8) the suspension viscosity n: 



„=|( 7 |AA- 1 | 7 ) 



[9] 



Eq.(9) is a key result, as it shows that all the singularities that 
appear near jamming within ASM are contained in a single 
operator, A/", which connects rheology to geometry. 



Analogy between elasticity and rheology 

The operator Af bears similarities with the well-studied stiff- 
ness matrix M that plays a central role in elasticity. By def- 
inition, the quadratic expansion of the energy in any elastic 
system is 5E — (SR\A4\SR) /2. Consider the case of an elas- 
tic network of un-stretched springs of unit stiffness. Labeling 
a spring by a, the energy can be expressed in terms of the 
spring elongation 8r a , namely: 

5E = J2 \ Sr « 2 = \{Sr\Sr) = ^{S8R\SSR) = ^(8R\S*S\8R) 

a 

[10] 

Together with the definition of Ai, Eq.(10) implies that 
M — 5*5. This expression is closely related to Af = 55*, 
building a connection between rheological properties of sus- 
pended hard particles and the elasticity of spring networks 
with identical geometry. In particular, the spectra of AA and 
Af are identical for positive eigenvalues A = cj 2 , where lo is 
referred to as the mode frequency. To see this let us assume 
that \8r(u})} is a normalized eigenvector of Af, then 

55*|<5r(w)) =Lj 2 \Sr(co)} . [11] 

Operating on the above equation with 5* we find: 

5'55'|<5r(w)) = A45*|5r(w)) = w 2 5*|<5r(o;)) , [12] 
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Fig. 3. The coordination displays scaling near jamming Sz = z c — z ~ in 
(a) two and (b) three dimensions. Particles barely connected to the contact network 
(making less than two contacts with their surrounding), the so-called rattlers, are 
removed from the analysis. Colors label different packing fractions (see Methods Sec- 
tion), growing from blue to red, and each data point pertains to a single configuration. 



implying that S t \Sr(u>)} is an eigenvector of M with the same 
eigenvalue uj 2 . Thus, there is a one-to-one correspondence 
between eigenvectors of M and N for positive eigenvalues. 

This observation suggests that the analytical techniques 
that were developed to study elasticity in random media ap- 
ply as well to rheology of suspensions, at least within ASM. 
Here we shall take another route by studying empirically the 
spectrum of N , identifying what the spectral signatures of the 
singularities are in the rheological laws, and compare those 
with the known singularity in the spectrum of M near the 
unjamming transition. 



Microscopic analysis of flow configurations 

A key variable that characterizes steric hindrance is the av- 
erage number of contacts per particle, or the coordination 
z = 2N C /N. This intuition is substantiated by Eq.(9), which 
indicates that the viscosity diverges when det(A/") = 0. Since 
J\f — SS 1 , J\f must display zero-modes if 5* does, which must 
occur if N c > ND, as 5* is of dimension ND x N c . For the 
coordination this argument implies that jamming must occur 
if z > 2D = z c , a necessary condition for rigidity derived 
by Maxwell [23]. Earlier numerical work [32] considered the 
scaling of the coordination with <f>, where it was found that 
for (f> < (f> c , Zo — z ~ $> c — <(> with zo = 3.8. However, this 
analysis is hampered by the rather large fluctuations of (j> c in 
finite size systems [35], in addition to errors originating from 
the inclusion of rattlers. Our numerical procedure allows for a 
careful measurement of the coordination, displayed in Fig. (3). 
We find instead: 



Sz 



1/P , 



[131 



with z c = 2D, S ~ 0.38 and 8 « 0.34 in two and three dimen- 
sions respectively. 



frequency spectra, or densities of states D{ui) defined as the 
number of modes per unit frequency per particle, are shown in 
Fig. (4) as jamming is approached. They consist of two struc- 
tures: a plateau of modes appearing above some frequency 
threshold ui* , and a gap at lower frequency, containing only 
one isolated mode of frequency o; m i n . Strikingly, this plateau 
is also present in the vibrational spectrum of simple amor- 
phous solids [14, 17], where its onset frequency scales asz-z c 
[15, 16]. To determine if this scaling law holds in configura- 
tions visited in flows as well, we rescale the frequency axis by 
Sz ~ l/p S - Fig. (4) displays a striking collapse of the plateau 
in the density of states, emphasizing the strong connection 
between elasticity near unjamming and the rheology of dense 
flows. 

The sole exception is the isolated mode at low-frequency 
that does not collapse with this rescaling. Rather, its fre- 
quency scales as aj m i n ~ l/p'~ with e « 0.51 independently of 
the spatial dimension, as shown in Fig. (5). Such bi-scaling is 
a new feature of flow not present in the elasticity of solids near 
unjamming, and it is of critical importance for the rheology. 

Using Eq. (9) we can separate the contributions of the first 
mode (Jo and of the plateau a — ao to the shear stress a. De- 
noting the lowest-frequency mode as \Sro), and the rest of the 
modes by \Sr u ), we obtain: 



£o7 <7l^o} 2 . fo7 



n 



+ 



E 



[141 



We find that the first term of Eq.(14), denoted <to, domi- 
nates stress near jamming: a — \ ctq. This stems from the 
fact that the first eigenvector has a very strong projection 
on a shear, that does not vanish as the number of particle 
increases: (7|<5ro)/(| Yl\ I x l l^ r o||) ~ 0.15 in three dimensions. 
Using that 1 1 — y 1 1 ~ N and that ||<5ro|| = 1 leads to rro ~ lA^mm- 
Since <to~o"~P we obtain w m in~l/^/p or e = 1/2, very close 
to the observation e~0.51. 

Note that the strong coupling between |<5ro) and I7) im- 
plies, together with Eq.(8), that |/}/||/|| — > \8ro) as p — > 00. 
Thus Fig. (la) and the movie in S.I., that display respectively 
an example of |/} and its evolution under shear very close to 
threshold, are very accurate representation of \5ro). 

The relative contribution of the rest of the density of 
states, corresponding to the second term in Eq.(14), vanishes 
as (<t — <To)/rj ~ p -0 ' 65 , as shown in Fig. (5). This exponent 
can also be rationalized by a simple scaling estimate, similar to 
arguments previously introduced for isotropic configurations 
near the unjamming transition [22, 33]. Assuming that the 
modes \Sru,) forming the plateau in D(to) have random orien- 
tations with respect to a shear [22] leads to {•y\8r UJ ) 2 ~ O(l). 
Thus 

I V (7l^) 2 1 v <«5rH| 7 ) 2 



era 



£07 



D(u)du) 

^2 



1/w* ~ l/Sz 



implying (a — ao)/a ~ Sz/a ~ 8z/p ~ p S1 which is in- 
deed p -0 ' 65 in three dimensions. Thus most of the spectrum 
contributes less and less to the divergence of stress as jam- 
ming is approached. This contribution may nevertheless play 
an important role in the corrections to scaling, such as those 
leading to a varying friction law in Eq.(2). We shall explore 
this possibility in future work. 



Density of states 

The symmetric operator M is readily obtained from the con- 
tact network, and its spectrum is computed numerically. The 



Discussion 

We have shown that the rheological properties of ASM flows 
are described by a single operator N , which is closely con- 
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Fig. 4. Spectral analysis of the operator governing flow as jamming is approached. The density of states D(id) averaged over 500 configurations for each pressure 
indicated in the legend, in (a) two and (b) three dimensions. The amplitude of the peak at low-frequency was rescaled to make it visible. Since it consists of one mode only, 
this amplitude vanishes in the thermodynamic limit. If the frequency axis is rescaled by the excess coordination p^ ~ Sz, as in (c,d), the emergence of a plateau in D(oj) 
collapses, indicating that the frequency at which this plateau appears follows UJ* ~ Sz. 




Fig. 5. Scaling of the lowest frequency UJ m i n vs. normalized pressure p in (a) 2 an 
by plotting the relative contribution (fT — <Tq)/<x of all the other modes to the stress 

nected to the stiffness matrix of elastic networks. This re- 
sult allows for a novel characterization of flow in terms of the 
spectrum of a single operator Af '. This spectrum presents re- 
markable features: it displays the plateau of modes controlling 
the anomalous elastic properties of amorphous solids near the 
unjamming transition, but also one mode at low-frequency 
responsible for the sharp divergence of the viscosity. Future 
work necessary to build a description of flow near jamming 
should explain (i) what configurations can generate such a bi- 
scaling spectrum (ii) why such configurations are selected by 
the dynamic and (iii) what fixes the relationship between cf> 
and z. 



d (b) 3 dimensions. This modes dominates the rheology near jamming, as can be seen 
that vanishes asp - 0-65. y^e relative contribution is found to be independent of N. 

It is likely that such a description will be of mean-field 
character. Near the unjamming transition the fact that a fre- 
quency scale uj* scales linearly with the excess coordination Sz 
in any dimension reflects the fact that spatial fluctuations of 
coordination are weak, and irrelevant. The unjamming tran- 
sition is in some sense a mean-field version of rigidity perco- 
lation where bonds are deposited randomly on a lattice and 
where fluctuations are important [34] , as a mean-field descrip- 
tion of the latter gives the correct elastic behavior near un- 
jamming [16]. Our finding that some frequency scale uj* also 
scales like Sz in ASM flows suggests that spatial fluctuations 
of coordination are irrelevant in dense suspensions too. This is 
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also supported by our observation that more generally, expo- 
nents appear to depend weakly, if at all, on spatial dimension. 

To conclude, we discuss the generality of our results. Con- 
stitutive relations in ASM appear quantitatively similar to 
experiments in real suspensions. This result might seem sur- 
prising at first glance, since the fluid is certainly strongly dis- 
turbed by flow near jamming, unlike what ASM assumes. We 
speculate why these two problems may fall in the same uni- 
versality class. In real flow dissipation is dominated by the 
viscous friction resulting from lubrication between neighbor- 
ing particles. The corresponding tangential friction coefficient 
has a weak (logarithmic) dependence on the gap between par- 
ticles, whereas the normal friction coefficient diverges as the 
inverse of the gap. However evidences [37, 27] support that 
due to the finite elasticity and/or roughness of the particles, 
real contacts are eventually formed [37, 27] . If so, one expects 
that the divergence of the normal friction coefficient will have 
a cut-off, and that the viscosity will be proportional to the 
square of the relative velocities between particles. In ASM, 
the viscosity is proportional to the square of the non-affinc 
velocities. However, we have checked that both quantities 
are approximatively proportional to each other as jamming is 
approached, supporting that this model is appropriate to de- 
scribe the scaling properties of suspension flows near jamming. 
Further support comes from a recent work on the viscoelas- 
ticity of amorphous solids near unjamming [33], showing that 
variations in the dissipation mechanism need not alter the 
scaling relations for the viscoelastic properties. 
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Methods 

The simulations are based on the ASM equations of motion, 
which prescribes a velocity to each particle for a given con- 
figuration. We propagate the system forward in time over 
small time steps, while carefully monitoring the formation of 
new contacts, or the destruction of existing contacts. For a 
complete description of the simulation employed see [31]. 

We simulated systems of JV2D = 4096 particles in two di- 
mensions, at packing fractions r/>2D = 0.82, 0.825, 0.83, 0.835, 
0.837, 0.838, 0.839, and 0.840. We also simulated systems 
of N30 = 1000 and 2000 particles in three dimensions, at 
packing fractions 3 d = 0.61, 0.625, 0.63, 0.635, 0.640, 0.642, 
and 0.643. In all simulations we used a square box with 
Lees-Edwards periodic boundary conditions [36] that allows 
for homogeneous shear flow profiles throughout the simula- 
tion cell. Our systems consist of a 50:50 mixture of small 
and large particles, where the ratio of the diameters d of 
the large and small particles is di a rge/d S maii = 1-4. We add 
a slight poly-dispersity of 3% in the particle sizes to avoid 
hexagonal patches in the two-dimensional systems [31]. Nor- 
mal modes are calculated with the LAPACK software package 
(http://www.netlib.org/lapack/). 
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